#########################
# Discounting Extreme Positions
# Analysis and replication file
# Laia Balcells, Sergi Martinez, and Ethan vanderWilden (corresponding author)
# ethan.vanderwilden@wisc.edu

# Last replicated: November 26, 2024
#########################

############ 1. Load in packages and data ####
# List of packages
pkg = c("tidyverse", "gridExtra", "haven", "modelsummary", "mediation", "Hmisc", 
        "collapse", "datawizard", "patchwork", "kableExtra", "vtable")

# Checks if they are installed, install if not
if (length(setdiff(pkg, rownames(installed.packages()))) > 0) {
  install.packages(setdiff(pkg, rownames(installed.packages())))
}
# Load
suppressMessages(suppressWarnings(lapply(pkg, library, character.only = TRUE)))
rm(list=ls())

setwd("C:/Users/19785/Dropbox/Ethan Sergi Laia/REPLICATION")
data <- read_sav("DiscountingData.sav")


############ 2. Write 'Helper' Functions  ####
# A. Get graphable coefficients for interaction model
# INPUT: Function takes in a dataset, the relevant outcome/dependent variable,
# the extreme position (LGBTQ or history), and whether control variables should be included in the model
# OUTPUT: Function outputs a dataframe with main coefficients and standard errors (easy to graph)
get_coefs <- function(data, DV = "sympathy", issue, controls = F){
  
  # If history issue, get regression, add results to the formatted output dataframe
  if(issue == "History"){
    if(controls == F){ myreg <- coef(summary(lm(scale(get(DV)) ~ normal + history + history*normal,  data = subset(data, data$treat %in% c(0:3)))))} 
    else {myreg <- coef(summary(lm(scale(get(DV)) ~ normal + history + history*normal + 
                                     age + sex + education + working + income + household_residents + pol_know + factor(region),
                                   data = subset(data, data$treat %in% c(0:3)))))}
    
    results <- data.frame(
      outcomes = c("Normalized", "Pro-Franco", "Normalized*\nPro-Franco"),
      estimate = myreg[c("normal", "history", "normal:history"), 1],
      se = myreg[c("normal", "history", "normal:history"), 2])
  } else{ #If LGBTQ+ issue position
    if(controls == F){ myreg <- coef(summary(lm(scale(get(DV)) ~ normal + lgbtq + lgbtq*normal,data = subset(data, data$treat %in% c(0:1, 4:5)))))} 
    else {myreg <- coef(summary(lm(scale(get(DV)) ~ normal + lgbtq + lgbtq*normal + 
                                     age + sex + education + working + income + household_residents + pol_know + factor(region),
                                   data = subset(data, data$treat %in% c(0:1,4:5)))))}
    
    results <- data.frame(
      outcomes = c("Normalized", "Anti-LGBTQ+", "Normalized*\nAnti-LGBTQ+"),
      estimate = myreg[c("normal", "lgbtq", "normal:lgbtq"), 1],
      se = myreg[c("normal", "lgbtq", "normal:lgbtq"), 2])
  }
  
  return(results)
}

# B. Get coefficient plot for main regression models
# INPUT: graphable coefficients (from get_coefs()), plus additional  graph specs
# OUTPUT: ggplot object with coefficient plot
coef_plot <- function(data, plot_title, axis_title = "\nStandardized effect", range = 0.5, jump = 0.25){
  
  #Ensure order of comparison stays consistent
  data$outcomes <- factor(data$outcomes, levels = data$outcomes)
  
  #plot 90 and 95% confidence intervals and coefficient point estimates
  myplot <- ggplot(data = data) +
    aes(x = estimate, y = outcomes) +
    geom_vline(xintercept = 0, size = 0.5)+
    geom_pointrange(aes(xmin= estimate-1.96*se, xmax= estimate+1.96*se), 
                    shape = 21, size = 0.8, linewidth = 1.8) +
    geom_pointrange(aes(xmin= estimate-1.645*se, xmax= estimate+1.645*se), 
                    shape = 21, size = 1.5, linewidth = 4, fill = "grey") +
    labs(x = axis_title,
         title = plot_title) +
    theme_bw()+
    theme(axis.text.y = element_text(size = 14, face = 'bold'),
          axis.text.x = element_text(size = 14),
          axis.title = element_text(size = 14, face = 'bold'),
          plot.title = element_text(size = 18, face = 'bold', hjust = 0.5),
          axis.title.y = element_blank(),
          panel.grid.minor = element_blank()) +
    scale_x_continuous(limits = c(-range, range), breaks = seq(-range, range, jump))+
    scale_y_discrete(limits=rev)
  
  #add line for manipulation check plot
  if(nrow(data)>3){ myplot <- myplot + geom_hline(yintercept = 2.5, linetype = "dashed")}
  
  return(myplot)
}


# C. Heterogeneous Splits
# INPUT: data subsets (ex: moderates and extremists), name of each dataset,
# outcome/dependent variable, extreme position (LGBTQ or history), and additional graph specs
# OUTPUT: ggplot object with coefficient plot separating results by heterogeneous sub-sample
het_splits <- function(data1, data2, name1, name2, dv = "sympathy", issue, plot_title, 
                       axis_title = "\nStandarized effect on sympathy for Vox", range = 0.5, jump = 0.25){
  
  #Get results for each subset-group
  results <- rbind(cbind(get_coefs(data1, DV = dv, issue = issue), group = rep(name1, 3)),
                  cbind(get_coefs(data2, DV = dv, issue = issue), group = rep(name2, 3)))
  
  ### Make plot
  results$outcomes <- factor(results$outcomes, levels = results$outcomes[1:3])
  results$group <- factor(results$group, levels = results$group[c(1, 4)])

  #plot 90 and 95% confidence intervals and coefficient point estimates
  myplot <- ggplot(data = results, aes(x = estimate, y = outcomes, color = group, fill = group)) +
    geom_vline(xintercept = 0, size = 0.5)+
  
    geom_pointrange(aes(xmin= estimate-1.96*se, xmax= estimate+1.96*se), 
                    shape = 21, size = 0.8, linewidth = 1.8,
                    position = position_dodge(width = -1/2)) +
    geom_pointrange(aes(xmin= estimate-1.645*se, xmax= estimate+1.645*se), 
                    shape = 21, size = 1.1, linewidth = 3,
                    position = position_dodge(width = -1/2)) +
    
    scale_color_manual(values = c("black", "darkgrey")) + scale_fill_manual(values = c('black', 'darkgrey')) +
    labs(x = axis_title, title = plot_title) +
    theme_bw()+
    theme(axis.text.y = element_text(size = 14, face = 'bold'),
          axis.text.x = element_text(size = 14),
          axis.title = element_text(size = 14, face = 'bold'),
          plot.title = element_text(size = 18, face = 'bold', hjust = 0.5),
          axis.title.y = element_blank(), legend.title = element_blank(),
          legend.text = element_text(size = 14),
          legend.position = c(0.85, 0.9), legend.box.background = element_rect(size = 1.1)) +
    scale_x_continuous(limits = c(-range, range), breaks = seq(-round(range, 1), round(range, 1), jump))+
    scale_y_discrete(limits=rev)
  
  return(myplot)
  
}


# D. Get coefficient plot for multiple outcomes (Sympathy, Ever Vote, Voting Vox) plot
# INPUT: dataset, extreme position (LGBTQ or history), graph spects
# OUTPUT: ggplot object with coefficient plot for multiple outcomes related to Vox
triple_outcomes <- function(data, issue, plot_title, axis_title = "\nStandarized effect on (outcome denoted in legend)", 
                            range = 0.5, jump = 0.25){
  
  #Get results
  results <- rbind( get_coefs(data, DV = "sympathy", issue = issue),  get_coefs(data, DV = "ever_vote", issue = issue),
    get_coefs(data, DV = "voteVox", issue = issue) )
  
  results$dv <- c(rep("Sympathy \n for Vox (0:100)", 3),  rep("Ever vote \n for Vox (1:5)", 3),  rep("Voting \n for Vox (0/1)", 3))
  
  
  ### Make plot
  results$outcomes <- factor(results$outcomes, levels = results$outcomes[1:3])
  results$dv <- factor(results$dv, levels = results$dv[c(1, 4, 7)])
  
  #plot 90 and 95% confidence intervals and coefficient point estimates
  myplot <- ggplot(data = results, aes(x = estimate, y = outcomes, color = dv, fill = dv)) +
    geom_vline(xintercept = 0, size = 0.5)+
    
    
    geom_pointrange(aes(xmin= estimate-1.96*se, xmax= estimate+1.96*se), 
                    shape = 21, size = 0.8, linewidth = 1.8,
                    position = position_dodge(width = -1/2)) +
    geom_pointrange(aes(xmin= estimate-1.645*se, xmax= estimate+1.645*se), 
                    shape = 21, size = 1.2, linewidth = 3,
                    position = position_dodge(width = -1/2)) +
    
    scale_color_manual(values = c("black", "darkgrey", "lightgrey")) +
    scale_fill_manual(values = c("black", "darkgrey", "lightgrey")) +
    
    labs(x = axis_title,
         title = plot_title) +
    theme_bw()+
    theme(axis.text.y = element_text(size = 14, face = 'bold'),
          axis.text.x = element_text(size = 14),
          axis.title = element_text(size = 14, face = 'bold'),
          plot.title = element_text(size = 16, face = 'bold', hjust = 0.5),
          axis.title.y = element_blank(), legend.title = element_blank(),
          legend.text = element_text(size = 14),
          legend.position = c(0.85, 0.8), legend.box.background = element_rect(size = 1.1)) +
    scale_x_continuous(limits = c(-range, range), breaks = seq(-range, range, jump))+
    scale_y_discrete(limits=rev)
  
  return(myplot)
  
}

############ 3. Main Results ####

### Figure 1 - LGBTQ results (render at 1400 x 500)
coef_plot(get_coefs(data, issue = "LGBTQ"), 
          "A. Sympathy for Vox", "\nStandardized effect on sympathy for Vox") +
coef_plot(get_coefs(data, DV = "M2.VoxHostile", issue = "LGBTQ"),
          "B. Vox is hostile towards LGBTQ+ people", "\nStandardized effect on perceptions of Vox authenticity") + 
  theme(axis.text.y = element_blank(), axis.ticks.y = element_blank())


### Figure 2 - Pro-Franco statement results (render at 1400 x 500)
coef_plot(get_coefs(data, issue = "History"), 
          "A. Sympathy for Vox", "\nStandardized effect on sympathy for Vox") +
coef_plot(get_coefs(data, DV = "M1.VoxReinstall", issue = "History"), 
          "B. Vox wants to re-install dictatorship", "\nStandardized effect on perceptions of Vox authenticity")+ 
  theme(axis.text.y = element_blank(), axis.ticks.y = element_blank())


### Figure 3 - differences between 'moderates' and 'extremists' (render at 1600 x 500)
het_splits(subset(data, data$ideology < 8), subset(data, data$ideology > 7), "Moderates (5:7)", "Far Right (8:10)",
           issue = "LGBTQ", plot_title = "A. Priming Anti-LGBTQ Statements") + 
het_splits(subset(data, data$ideology < 8), subset(data, data$ideology > 7), "Moderates (5:7)", "Far Right (8:10)",
           issue = "History", plot_title = "B. Priming Pro-Franco Statements")


############## 4. Appendix ####

#### A.1.1 Demographic representativeness ####

#Because demographic categories are categorical, we have to convert data into binaries
repdata <- data #save data as new object

#Get percentages for each variable
repdata$M1 <- data$sex
repdata$AG1 <- ifelse(repdata$age %in% c(18:25), 1, 0)
repdata$AG2 <- ifelse(repdata$age %in% c(26:35), 1, 0)
repdata$AG3 <- ifelse(repdata$age %in% c(36:45), 1, 0)
repdata$AG4 <- ifelse(repdata$age %in% c(46:55), 1, 0)
repdata$AG5 <- ifelse(repdata$age %in% c(56:65), 1, 0)
repdata$AG6 <- ifelse(repdata$age >65, 1, 0)
repdata$R1 <- ifelse(repdata$region == 1, 1, 0) #Andalucia
repdata$R2 <- ifelse(repdata$region == 2, 1, 0) #Aragon
repdata$R3 <- ifelse(repdata$region == 3, 1, 0) #Asturias
repdata$R4 <- ifelse(repdata$region == 4, 1, 0) #Islas Baleares
repdata$R5 <- ifelse(repdata$region == 5, 1, 0) #Islas Canarias
repdata$R6 <- ifelse(repdata$region == 6, 1, 0) #Cantabria
repdata$R7 <- ifelse(repdata$region == 8, 1, 0) #note: swapped in INE data, so make Castilla y Leon 8th
repdata$R8 <- ifelse(repdata$region == 7, 1, 0) #Castilla la Mancha
repdata$R9 <- ifelse(repdata$region == 9, 1, 0) #Catalunya
repdata$R10 <- ifelse(repdata$region == 10, 1, 0) #Valencia
repdata$R11 <- ifelse(repdata$region == 11, 1, 0) #Extremadura is 11
repdata$R12 <- ifelse(repdata$region == 12, 1, 0) #Galicia
repdata$R13 <- ifelse(repdata$region == 13, 1, 0) #Madrid
repdata$R14 <- ifelse(repdata$region == 14, 1, 0) #Murcia
repdata$R15 <- ifelse(repdata$region == 15, 1, 0) #Navarra
repdata$R16 <- ifelse(repdata$region == 16, 1, 0) #Pais Vasco
repdata$R17 <- ifelse(repdata$region == 17, 1, 0) #La Rioja

#Set up dataframe to store INE statistics + sample statistics
#On gender; age group; and region
rep.tab <- data.frame(
  row.names = c("Male", "18-25", "26-35", "36-45", "46-55", "56-65", "66+",
                "Andalucía", "Aragón", "Asturias", "las Islas Baleares", "las Islas Canarias",
                "Cantabria", "Castilla y León", "Castilla-La Mancha", "Catalunya", "Valencia",
                "Extremadura", "Galicia", "Madrid", "Murcia", "Navarra", "País Vasco", "La Rioja"),
  
  sample.prop = rep(NA, 24), 
  census.prop = c(0.490, 0.097, 0.131, 0.169, 0.187, 0.158, 0.256,
                  0.179, 0.028, 0.021, 0.025, 0.046, 0.012, 0.050, 0.043,
                  0.165, 0.109, 0.022, 0.056, 0.143, 0.032, 0.014, 0.046, 0.007),
  difference = rep(NA, 24),
  p.value = rep(NA, 24)
)

#For eadh variable, gather sample proportion, difference from INE; and p-value for t.test
for (i in 1:24){
  rep.tab$sample.prop[i] <- round(mean(repdata[, i+33][[1]], na.rm = T), 3)
  rep.tab$difference[i] <- round(abs(rep.tab$sample.prop[i] - rep.tab$census.prop[i]), 3)
  rep.tab$p.value[i] <- round(t.test(repdata[,i+33], mu = rep.tab$census.prop[i], alternative = 'two.sided')$p.value, 3)
}

#Print table (latex)
kbl(rep.tab, format = 'latex', booktabs = T)






#### A.1.1 Comparing to sample demographics to CIS ####

#Load in CIS data, September 2023 (3420)
MD3420 <- read_sav("additional_data/cis2023sept3420.sav")

#Clean data and make numeric [codebook in replication files]
MD3420$educacion <- as.numeric(MD3420$NIVELESTENTREV)
MD3420$educacion <- ifelse(MD3420$educacion==98, NA, MD3420$educacion)
MD3420$sexo <- (as.numeric(MD3420$SEXO)-1)
MD3420$ingresos <- as.numeric(MD3420$INGRESHOG)
MD3420$working <- ifelse(MD3420$SITLAB==1 | MD3420$SITLAB==7 | MD3420$SITLAB==8, 1, 0)

## TABLE A2 - full sample
suppressWarnings(sumtable(MD3420, vars = c("sexo", "EDAD", "educacion", "working"),
                          labels = c("Gender", "Age", "Education", "Working"), out='latex'))

## TABLE A3 - right-only subsample
suppressWarnings(sumtable(MD3420[MD3420$ESCIDEOL>4 & MD3420$ESCIDEOL<12,], 
                         vars = c("sexo", "EDAD", "educacion", "working"),
                         labels = c("Gender", "Age", "Education", "Working"), out='latex'))


#### A.1.3 Descriptive statistics Table ####


#Store variable names needed for descriptive statistics
variables <- c('sympathy', 'ever_vote', 'voteVox', 'M1.VoxReinstall', 'M2.VoxHostile',
               'vote_guess', 'coalition.manip',
               'sex', 'age', 'education', 'household_residents', 'income', 'working', 
               'SP_nationalism', 'CCAA_nationalism', 'ideology', 'pol_know')

# TABLE A4 (latex)
sumtable(data, vars = variables, labels = variables,
         numformat = formatfunc(digits = 2, nsmall = 2, big.mark = ','), out = 'latex')



#### A.1.3 Outcome distributions ####

## FIGURE A2 (a, b, c) - Vox affinities
#Get means for graph labels
mean(data$sympathy, na.rm = T) #43.01
mean(data$ever_vote, na.rm = T) #2.53
mean(data$voteVox, na.rm = T) #0.178

#Vox support
sym <- ggplot(data = data) +
  geom_histogram(aes(x = sympathy), color = "black", fill = "darkgreen", alpha = 0.7, bins = 25) +
  theme_bw() + geom_vline(xintercept = 43.01, color = "black") +
  labs(title = "Sympathy for Vox (Mean response: 43.01)", x = "Sympathy score", y = "Frequency") +
  scale_x_continuous(limits = c(-5, 105), breaks = seq(0, 100, 10)) +
  theme(text = element_text(size = 12), plot.title = element_text(size = 16, face = 'bold', hjust = 0.5),
        panel.grid = element_blank())

ever <- ggplot(data = data) +
  geom_bar(aes(x = ever_vote), color = "black", fill = "darkgreen", alpha = 0.7) +
  theme_bw() + geom_vline(xintercept = 2.53, color = "black") +
  labs(title = "Consider ever voting Vox (Mean response: 2.53)", 
       x = "1 = Never vote Vox, 5 = Will definitely vote Vox in the future", y = "Frequency") +
  scale_x_continuous(limits = c(0.49, 5.51), breaks = seq(1, 5, 1)) +
  theme(text = element_text(size = 12), plot.title = element_text(size = 16, face = 'bold', hjust = 0.5),
        panel.grid = element_blank())

vote <- ggplot(data = data) +
  geom_bar(aes(x = voteVox), color = "black", fill = "darkgreen", alpha = 0.7) +
  theme_bw() + labs(title = "Planning to vote Vox (Mean response: 0.178)", 
                    x = "0 = Other party/no Vote, 1 = Voting Vox", y = "Frequency") +
  scale_x_continuous(limits = c(-0.49, 1.51), breaks = seq(0, 1, 1)) +
  theme(text = element_text(size = 12), plot.title = element_text(size = 16, face = 'bold', hjust = 0.5),
        panel.grid = element_blank())

#render figure A2 (1650 x 450)
sym + ever + vote

## FIGURE A3 (a, b) - authenticity
#Get means for graph labels
mean(data$M2.VoxHostile, na.rm = T) # 3.28
mean(data$M1.VoxReinstall, na.rm = T) # 2.38

#Vox is hoslite towards LGBTQ+ community
hostileLBGTQ <- ggplot(data = data) +
  geom_bar(aes(x = M2.VoxHostile), color = "black", fill = "darkgreen", alpha = 0.7) +
  theme_bw() + geom_vline(xintercept = 3.28, color = "black") +
  labs(title = "Vox is hostile to LGBTQ+ people (Mean response: 3.28)", 
       x = "1 = Strongly disagree, 5 = Strongly agree", y = "Frequency") +
  scale_x_continuous(limits = c(0.49, 5.51), breaks = seq(1, 5, 1)) +
  theme(text = element_text(size = 12), plot.title = element_text(size = 16, face = 'bold', hjust = 0.5),
        panel.grid = element_blank())

#Vox would re-install dictatorship
reinstall <- ggplot(data = data) +
  geom_bar(aes(x = M1.VoxReinstall), color = "black", fill = "darkgreen", alpha = 0.7) +
  theme_bw() + geom_vline(xintercept = 2.38, color = "black") +
  labs(title = "Vox would re-install dictatorship (Mean response: 2.38)",
       x = "1 = Strongly disagree, 5 = Strongly agree", y = "Frequency") +
  scale_x_continuous(limits = c(0.49, 5.51), breaks = seq(1, 5, 1)) +
  theme(text = element_text(size = 12), plot.title = element_text(size = 16, face = 'bold', hjust = 0.5),
        panel.grid = element_blank())

#Render figure A3 (1250 x 450)
hostileLBGTQ + reinstall



#### A.1.4 Balance Tables ####

#Save variables (pre-treatment) for balance checks
variables <- c('sex', 'age', 'education', 'household_residents', 'income', 'working', 
               'SP_nationalism', 'CCAA_nationalism', 'ideology', 'pol_know')

#Create dataframe to store results
results <- data.frame("Variable" = variables, "Low" = rep(NA, length(variables)),
                      "High" = rep(NA, length(variables)), "SC" = rep(NA, length(variables)),
                      "NC" = rep(NA, length(variables)), "SH" = rep(NA, length(variables)),
                      "NH" = rep(NA, length(variables)), "SL" = rep(NA, length(variables)),
                      "NL" = rep(NA, length(variables)), "ANOVA" = rep(NA, length(variables)))

#for each variable, record range, means, and anova p-value for balance across treatment arms
for (i in 1:length(variables)){
  info <- c()
  info <- append(info, range(data[, c(variables[i])][[1]], na.rm = T)) # range
  for(j in 0:5){info <- append(info, mean(subset(data, treat == j)[ , c(variables[i])][[1]], na.rm = T))} #means
  info <- append(info, as.numeric(anova(lm(get(variables[i]) ~ factor(treat), data = data))[1,5])) # ANOVA p-value
  
  results[i, 2:10] <- info
}

#Round numeric results
results[,2:10] <- round(results[, 2:10], 3)

#TABLE A5 - Print results
kbl(results, format = 'latex', booktabs = T)



#### A.2.1 Manipulation checks ####

#create dataframe to store results
results <- data.frame( outcomes = c("Perceived accept-\nability of coalitions\n(effect of 'normalized')", 
                                    "Vote approx. for Vox\n(effect of 'normalized')", 
                                    "Vox is hostile\ntowards LGBTQ+ people\n(effect of anti-LGBTQ+ statement)", 
                                    "Vox would re-install autocracy\n(effect of pro-Franco statement)"),
               estimate = rep(NA, 4), se = rep(NA, 4))

#A. Normalization (keep main effect of normalized)
results[1, 2:3] <- coef(summary(lm(scale(coalition.manip) ~ normal + lgbtq + history +
                                     lgbtq*normal + history*normal, data = data)))[2, 1:2]
results[2, 2:3] <- coef(summary(lm(scale(vote_guess) ~ normal + lgbtq + history + 
                                     lgbtq*normal + history*normal, data = data)))[2, 1:2]


#B. LGBT (keep main effects of LGBTQ statement)
results[3, 2:3] <- coef(summary(lm(scale(M2.VoxHostile) ~ lgbtq + normal + normal*lgbtq, 
                                   data = subset(data, data$treat %in% c(0:1, 4:5)))))[2, 1:2]

#C. Authoritarian history (keep main effect of History statement)
results[4, 2:3] <- coef(summary(lm(scale(M1.VoxReinstall) ~ history  + normal + history*normal, 
                                   data = subset(data, data$treat %in% c(0:3)))))[2, 1:2]

## Render Figure A4 (900 x 550)
coef_plot(results, "Manipulation Checks", "\nStandardized effect of relevant manipulation") +
  theme(axis.text.y = element_text(face = 'plain'))


#### A.2.2 Did party normalization normalize the issue positions? ####

#Check if normalization prime led to change on expression of issue support 
# NOTE: control group, no issue prime
mod1 <- summary(lm(M3.FrancoGood ~ normal, data = subset(data, data$treat %in% c(0:1))))
mod2 <- summary(lm(M4.DoMoreLGBTQ ~ normal, data = subset(data, data$treat %in% c(0:1))))

#TABLE A6 - print regression outputs (cleaned up in paper)
suppressWarnings(modelsummary(list(mod1, mod2), output = 'latex', stars = T))


#### A.2.3 Tabular results from main figures ####
#Make treatment assignment for Vox statement a factor variable
data$position <- factor(data$position, levels = c('null', 'lgbtq', 'history'))

#Record main regressions
main <- summary(lm(sympathy ~ position + normal + position*normal, data = data))
mL <-summary(lm(M2.VoxHostile ~ position + normal + position*normal, data = subset(data, position != 'history')))
mH <-summary(lm(M1.VoxReinstall ~ position + normal + position*normal, data = subset(data, position != 'lgbtq')))
mods <- summary(lm(sympathy ~ position + normal + position*normal, data = subset(data, ideology <8)))
exts <- summary(lm(sympathy ~ position + normal + position*normal, data = subset(data, ideology >7)))

#TABLE A7 - print regressions, cleaned up in paper
suppressWarnings(modelsummary(list(main, mL, mH, mods, exts), stars = T, output = 'latex', booktabs = F))


#### A.2.4 Examining more outcomes related to Vox ####

# Figure A5 render at 1600 x 500
triple_outcomes(data, issue = "LGBTQ", plot_title = "A. Normalization and anti-LGBTQ+ statements") +
triple_outcomes(data, issue = "History", plot_title = "B. Normalization and pro-Franco statements")


#### A.2.5 Additional analyses for Moderates and extremists ####

# Figure A6, Mechanisms (discounted authenticity): render at 1600 x 500
het_splits(subset(data, data$ideology < 8), subset(data, data$ideology > 7), "Moderates (5:7)", "Far Right (8:10)",
           dv = "M2.VoxHostile", issue = "LGBTQ", plot_title = "A. Priming Anti-LGBTQ+ Statements",
           axis_title = "\nStandarized effect on perceptions of Vox authenticity") +
het_splits(subset(data, data$ideology < 8), subset(data, data$ideology > 7), "Moderates (5:7)", "Far Right (8:10)",
           dv = "M1.VoxReinstall", issue = "History", plot_title = "B. Priming Pro-Franco Statements",
           axis_title = "\nStandarized effect on perceptions of Vox authenticity")


# Figure A7, Alternative cutoff points: render at 1600 x 500
het_splits(subset(data, data$ideology < 9), subset(data, data$ideology > 8), "Moderates (5:8)", "Far Right (9:10)",
           issue = "LGBTQ", plot_title = "A. Priming Anti-LGBTQ+ Statements", range = 0.6, jump = 0.3) +
het_splits(subset(data, data$ideology < 9), subset(data, data$ideology > 8), "Moderates (5:8)", "Far Right (9:10)",
                         issue = "History", plot_title = "B. Priming Pro-Franco Statements", range = 0.6, jump = 0.3)


# Figure A8, Vox and Non-vox voters: render at 1600 x 500
het_splits(subset(data, data$voteVox == 0), subset(data, data$voteVox == 1), "Non-Vox voters", "Vox voters",
           dv = "sympathy", issue = "LGBTQ", plot_title = "A. Priming Anti-LGBTQ+ Statements", range = 0.54) +
het_splits(subset(data, data$voteVox == 0), subset(data, data$voteVox == 1), "Non-Vox voters", "Vox voters",
                         dv = "sympathy", issue = "History", plot_title = "B. Priming Pro-Franco Statements", range = 0.54)

#### A.2.6 Mediation ####
# NOTE: not looking at interaction with mediation (gets severely underpowered)...
# just looking at subset treated with Vox statement
# Examine split with moderates and extreme right (most interesting findings from the results section... 
# driving later analysis choices [note: re-specifying with full subset, approaching sig with ACME])
set.seed(1234)


### More moderate respondents
med1 <- subset(data, data$treat %in% c(4:5) & data$ideology<8)
med1$sympathy <- scale(med1$sympathy)

### 1. Is there discounting for Vox's anti-LGBT statements?
#fit model where treatment predicts mediator
med.fit <- lm(M2.VoxHostile ~ normal , data = med1) # significant in a 1 tail test)

#fit model predicting outcome based on treatment and mediator
out.fit <- lm(sympathy ~ M2.VoxHostile + normal, data = med1)

#simulation to calculate average causal mediation effect
med.out <- mediate(med.fit, out.fit, boot = TRUE, treat = "normal", mediator = "M2.VoxHostile", sims = 1000)
summary(med.out)

#Add results to table
med_results <- cbind(as.data.frame(t(as.matrix(data.frame(
  ACME = c( 0.1165, 0.0124, 0.22) , ADE = c(0.1519, 0.0346, 0.27) ,  `Total Effect` = c(0.2684, 0.1094, 0.42),
  row.names = c("Estimate", "CI_low", "CI_high") )))), Coefs = c("ACME", "ADE", "Total Effect"),
  group = rep("Moderates", 3))

#test sensitivity of rho parameter
sens.out <- medsens(med.out, rho.by = 0.05, effect.type = "indirect",
                    sims = 1000)
summary(sens.out) #rho at which ACME goes to 0: -0.65

#plot sensitivity (FIGURE A9b)
plot(sens.out, sens.par = "rho", main = "Sensitivity for Discounting Mechanism (Moderate sub-sample)", 
     ylim = c(-0.5, 0.5), xlim = c(-0.9, 0.2))


### More Extreme Respondents
med1 <- subset(data, data$treat %in% c(4:5) & data$ideology>7)
med1$sympathy <- scale(med1$sympathy)

### 1. Is there discounting for Vox's anti-LGBT statements?
#fit model where treatment predicts mediator
med.fit <- lm(M2.VoxHostile ~ normal , data = med1) # significant in a 1 tail test)

#fit model predicting outcome based on treatment and mediator
out.fit <- lm(sympathy ~ M2.VoxHostile + normal, data = med1)

#simulation to calculate average causal mediation effect
med.out <- mediate(med.fit, out.fit, boot = TRUE, treat = "normal", mediator = "M2.VoxHostile", sims = 1000)
summary(med.out)

#Add results to table
med_results2 <- cbind(as.data.frame(t(as.matrix(data.frame(
  ACME = c(0.0293, -0.1054, 0.17) , ADE = c(0.0326, -0.1770, 0.25) ,  `Total Effect` = c(0.0619, -0.1928, 0.33),
  row.names = c("Estimate", "CI_low", "CI_high") )))), Coefs = c("ACME", "ADE", "Total Effect"),
  group = rep("Far Right", 3))

#test sensitivity of rho parameter
sens.out <- medsens(med.out, rho.by = 0.05, effect.type = "indirect",
                    sims = 1000) #rho at which ACME goes to 0: -0.55


#plot sensitivity (Not shown, bc no unobserved confounding needed to render insignificant)
# e.g., already insignificant ACME
plot(sens.out, sens.par = "rho", main = "Sensitivity for Discounting Mechanism (Far Right Respondents)", 
     ylim = c(-0.5, 0.5), xlim = c(-0.9, 0.2))


#bind mediation results for moderates and extremists (show on one plot)
med <- rbind(med_results, med_results2)


#Plot mediation coefficients (FIGURE A9a)
ggplot(data = med, aes(x = Estimate, y = Coefs, color = group)) +
  geom_pointrange(aes(xmin = CI_low, xmax = CI_high), 
                  size = 0.8, linewidth = 2,
                  position = position_dodge(width = -1/2)) +
  geom_vline(xintercept = 0, linetype = "dashed") +
  scale_color_manual(values = c("black", "darkgrey"))+
  
  labs(y = "", x = "Standardized estimate", title = "Mediation effects, discounting mechanism") +
  scale_y_discrete(limits = rev) + scale_x_continuous(limits = c(-0.2, 0.5), breaks = seq(-0.2, 0.4, 0.2)) +
  theme_bw() +
  theme(plot.title = element_text(size = 16),
        legend.text = element_text(size = 12),
        legend.title = element_blank(),
        legend.position = "bottom",
        axis.text = element_text(size = 14),
        axis.title.x = element_text(size = 14))



#### A.2.7 Heterogeneous effects, demographics ####

### Figure A10, LGBTQ (render at 2000 x 2000)
gender <- het_splits(subset(data, data$sex == 0),subset(data, data$sex == 1), 
                     "Women", "Men",issue = "LGBTQ", plot_title = "Gender")
age <- het_splits(subset(data, data$age>52.1),subset(data, data$age <51.1), 
                  "Over 51", "51 or younger", issue = "LGBTQ", plot_title = "Age") +
  theme(axis.text.y = element_blank(), axis.ticks.y = element_blank())
income <- het_splits(subset(data, data$income > 6.01 ), subset(data, data$income < 6.01), 
                     "Income > 2500 Euro/month", "Income < 2500 Euro/month", issue = "LGBTQ", plot_title = "Income")
sp_id <- het_splits(subset(data, data$SP_nationalism > 8.3 ), subset(data, data$income < 8.3), 
                    "High Sp.ID (9:10)", "Low Sp.ID (<9)", issue = "LGBTQ", plot_title = "Spanish Nationalism") +
  theme(axis.text.y = element_blank(), axis.ticks.y = element_blank())
ccaa_id <- het_splits(subset(data, data$CCAA_nationalism > 8.3 ),subset(data, data$CCAA_nationalism < 8.3),
                      "High CCAA.ID (9:10)", "Low CCAA.ID (<9)", issue = "LGBTQ", plot_title = "Regional Identification")
pol_k <- het_splits(subset(data, data$pol_know == 1),subset(data, data$pol_know != 1), 
                    "Correct", "Incorrect", issue = "LGBTQ", plot_title = "Political Knowledge") +
  theme(axis.text.y = element_blank(), axis.ticks.y = element_blank())


gender + age + income + sp_id + ccaa_id + pol_k + plot_layout(ncol = 2)



### Figure A11, HISTORY (render at 2000 x 2000)
gender <- het_splits(subset(data, data$sex == 0),subset(data, data$sex == 1), 
                     "Women", "Men", issue = "History", plot_title = "Gender")
age <- het_splits(subset(data, data$age>52.1),subset(data, data$age <51.1), 
                  "Over 51", "51 or younger", issue = "History", plot_title = "Age")+
  theme(axis.text.y = element_blank(), axis.ticks.y = element_blank())
income <- het_splits(subset(data, data$income > 6.01 ), subset(data, data$income < 6.01), 
                     "Income > 2500 Euro/month", "Income < 2500 Euro/month", issue = "History", plot_title = "Income")
sp_id <- het_splits(subset(data, data$SP_nationalism > 8.3 ),subset(data, data$income < 8.3), 
                    "High Sp.ID (9:10)", "Low Sp.ID (<9)", issue = "History", plot_title = "Spanish Nationalism")+
  theme(axis.text.y = element_blank(), axis.ticks.y = element_blank())
ccaa_id <- het_splits(subset(data, data$CCAA_nationalism > 8.3 ), subset(data, data$CCAA_nationalism < 8.3), 
                      "High CCAA.ID (9:10)", "Low CCAA.ID (<9)", issue = "History", plot_title = "Regional Identification")
pol_k <- het_splits(subset(data, data$pol_know == 1),subset(data, data$pol_know != 1), 
                    "Correct", "Incorrect", issue = "History", plot_title = "Political Knowledge")+
  theme(axis.text.y = element_blank(), axis.ticks.y = element_blank())

gender + age + income + sp_id + ccaa_id + pol_k + plot_layout(ncol = 2)




#### A.3 Vox ideological profile ####

#Read in manifesto project data
dat <- read_dta("additional_data/MPDataset_MPDS2023a_stata14.dta", encoding = 'UTF-8')

## FIGURE A12a (within Spain comparison)
dat <- dat[dat$countryname=="Spain",]
dat <- dat[dat$partyabbrev=="VOX" | dat$partyabbrev=="C's" | dat$partyabbrev=="PSOE" | dat$partyabbrev=="PP" | dat$partyabbrev=="IU" | dat$partyabbrev=="UP", ]

manifestoSpain <- dat[dat$countryname=="Spain" & dat$coderyear>2016,] %>%
  ggplot(aes(x = rile, y = per603, colour= partyabbrev))+
  geom_jitter(alpha=1, aes(size = pervote)) +
  geom_text(data = dat[dat$edate=="2019-11-10",], aes(label=partyabbrev), size=6, vjust = -.9, hjust = .7, family= "serif") +
  labs(x="Left-right scale",
       y="Traditional\nmorality") +
  scale_colour_manual(name = "Political parties:", 
                      values = c("orange", "purple", "lightblue", "red", "green")) +
  theme_bw(base_size = 17) +
  theme(text = element_text(family="serif"), panel.grid = element_blank(), 
        axis.title.y = element_text(angle = 0), legend.position = "none")


## FIGURE A12b (across far right party comparison)
#Re-read data for radical right subset
dat <- read_dta("additional_data/MPDataset_MPDS2023a_stata14.dta", encoding = 'UTF-8')

dat <- dat[dat$parfam==70,]
dat$year=substr(dat$edate,1,4)
dat$vox <- ifelse(dat$partyabbrev=="VOX", 1, 0)
dat$vox <- factor(dat$vox, levels = c(0,1), labels = c("Other far-right parties", "Vox"))

manifestoFarRight <- dat[dat$coderyear>2010,] %>%
  ggplot(aes(x = rile, y = per603, colour= vox))+
  geom_jitter(alpha=1, aes(size = pervote)) +
  geom_text(data = dat[dat$year>2016 & dat$year<2020,], aes(label=partyabbrev), size=6, vjust = -.9, hjust = .7, family= "serif") +
  labs(x="Left-right scale",
       y="Traditional\nmorality") +
  scale_colour_manual(values = c("gray", "black")) +
  theme_bw(base_size = 17) +
  theme(text = element_text(family="serif"), panel.grid = element_blank(), 
        axis.title.y = element_text(angle = 0), legend.position = "none")


## Figure A12 (render at 1350 x 600)
manifestoSpain + manifestoFarRight
